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Gaussian boson sampling (GBS) is not only a feasible protocol for demonstrating quantum computational 
advantage, but also mathematically associated with certain graph-related and quantum chemistry problems. In 
particular, it is proposed that the generated samples from the GBS could be harnessed to enhance the classical 
stochastic algorithms in searching some graph features. Here, we use Jilizhdng, a noisy intermediate-scale 
quantum computer, to solve graph problems. The samples are generated from a 144-mode fully-connected 
photonic processor, with photon-click up to 80 in the quantum computational advantage regime. We investigate 
the open question of whether the GBS enhancement over the classical stochastic algorithms persists—and how it 
scales—with an increasing system size on noisy quantum devices in the computationally interesting regime. We 
experimentally observe the presence of GBS enhancement with large photon-click number and a robustness of 
the enhancement under certain noise. Our work is a step toward testing real-world problems using the existing 
noisy intermediate-scale quantum computers, and hopes to stimulate the development of more efficient classical 


and quantum-inspired algorithms. 


Recent experiments have constructed noisy intermediate- 
scale quantum (NISQ) devices and shown increasingly more 
convincing evidence for quantum computational advantage 
OHS], a milestone that demonstrates that the quantum devices 
can solve sampling problems overwhelmingly faster than clas- 
sical computers. A natural next step is to test whether these 
NISQ devices can solve problems of practical interest. 

Gaussian boson sampling (GBS) [6], a variant of the origi- 
nal Aaronson-Arkipov boson sampling [7], has attracted con- 
siderable attention for its potential applications in graph- 
related problems, quantum chemistry, and machine learning 
[BHI4]. This is because the GBS has underlying mathemat- 
ical connection with graph theory. Therefore, theories [9 
suggested that the generated samples from the GBS might 
give enhancement over classical stochastic algorithms in solv- 
ing some graph problems. Moreover, the highly connected 
topology of the GBS photonic processor can naturally address 
problems on nonplanar graphs [15]. 

Proof-of-principle demonstrations of solving graph prob- 
lems assisted by the GBS have been reported [161/19], how- 
ever, in regimes where the GBS devices dynamics can be 
easily simulated on classical computers. An important and 
open question is whether the GBS could give enhancement 
on increasingly larger devices in the computationally interest- 
ing regime, and how the performance is affected by noise in 
NISQ devices. Furthermore, previous demonstrations on find- 
ing dense subgraphs could only address the problem with non- 
negative-valued sampling matrices, for which efficient clas- 


sical algorithms of estimating the sampling probability exist 
and a quantum-inspired classical algorithm was re- 
cently developed [22]. 

Here, we test solving nonplanar graph problems on the 
NISQ photonic quantum processor, Jilizhdng, with 50 single- 
mode squeezed states input into a 144-mode linear optical 
network [2] B]. We operate Jitizhang in the computationally 
interesting regime to enhance stochastic algorithms solving 
two graph problems, namely the Max-Haf problem [9] and 
the dense k-subgraph problem [10]. We benchmark how the 
performance scales as a function of the GBS size, and how it 
is influenced by certain noise [23]. 

In the GBS, arrays of squeezed vacuum states are sent 
through a multi-mode interferometer and sample the output 
scattering events. Due to its Gaussian properties, the output 
state can be described by its Husimi covariance matrix og 
(24) (25], for which the sampling matrix is expressed as 


JOZ —1 
f l (1-03'). (1) 
The sampling matrix A is in a block matrix form 
A L 
A= | | | 2) 


where A is a symmetric matrix, and L = 0 if the Gaussian 
state is a pure state. 


log,o{|Hafnian|] 


-13 -11 -9 -7 -5 -3 
log4o[Torontonian] 


l0g4o[Density] 


“43 mn 3 7 5 3 
log,,{Torontonian] 


FIG. 1. Principle of the experiment. (a) The correspondence between a GBS experiment and a graph. Each output port of the GBS corresponds 
to a vertex of the graph, and the clicked ports correspond to a subgraph whose vertices and edges are marked in red. (b) The Monte Carlo 
simulation results showing Hafnian’s dependence on Torontonian for the randomly sampled four-mode complex-valued matrix. (c) The Monte 
Carlo simulation results showing graph density’s dependence on Torontonian for the randomly sampled four-mode complex-valued matrix. 


An illustration of the correspondence between a graph and a 
GBS setup is shown in Fig. Ma. Any undirected graph can be 
represented by its adjacency matrix A which is symmetric— 
i.e. A;; = Aj;—and the adjacency matrix element A;; cor- 
responds to the weighted value of the edge connecting vertex 
i to vertex j. The adjacency matrix can be encoded into the 
sampling matrix A of a pure state GBS with a proper rescaling 
factor c (see the Supplemental Material [26]): 


cA 0 
wale 


and by Takagi-Autonne decomposition the corresponding 
GBS setup can be constructed. Each mode of the output light 
field maps to a column and row of the adjacency matrix, and 
each GBS sample corresponds to a submatrix of the sampling 
matrix A by taking the elements of the corresponding rows 
and columns. Once the relationship between the GBS device 
and the adjacency matrix of the graph under study is estab- 
lished, the GBS samples, whose probability is positively cor- 
related to the mathematical quantity called Torontonian 
(for threshold detection) or Hafnian (for photon-number- 
resolving detection) of the corresponding submatrix, are har- 
nessed to enhance solving the graph problems of interest. 

We study the GBS enhancement on solving the Max-Haf 


problem and dense k-subgraph problem. The Max-Haf prob- 
lem is, for a complex-valued matrix B of any dimension, to 
find a submatrix Bs of fixed even dimension k = 2m, with 
the largest Hafnian in square of absolute value. Hafnian was 
originally introduced in interacting quantum field theory and 
plays a variety of roles in physics and chemistry [8] [30}/36). 
When the matrix is an adjacency matrix composed of Os and 
1s, Hafnian can be interpreted as the number of perfect match- 
ing of the graph [11]. The Max-Haf problem is known to be- 
long to the NP-hard complexity class [9]. 

The dense k-subgraph problem is, for an n-vertex graph 
G with adjacency matrix A, to find its subgraph of k < n 
vertices G's with the largest density 


S (As);5|5 (4) 


where Ag is the adjacency matrix of G's. 

The dense k-subgraph problem is of fundamental interest 
in both mathematics and applied fields like data mining, 
bioinformatics, finance and network analysis [B8445]. Al- 
though there are deterministic algorithms for finding subgraph 
of large density, they can be fooled and thus stochastic algo- 
rithms are important in certain scenarios [10]. 
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FIG. 2. The GBS enhancement on Max-Haf and dense k-subgraph problem solving. (a,b) The GBS enhancement on finding the 12-vertex 
subgraph with the largest Hafnian out of a 144-vertex graph for SA (a) and RS (b) algorithm. The experiment obtains a mean photon-click 
number of 21. A set of 12-photon-click samples are post-selected for this study. The X-axis corresponds to the number of steps used in 
the optimization process, while the Y-axis indicates the corresponding largest Hafnian found. Each curve is the mean largest Hafnian out 
of 100 trials, and the shaded area indicates standard error. (c) The graph corresponding to the experiment, and the subgraph of maximum 
Hafnian found by the GBS-enhanced SA algorithm (vertices marked in red). 119 unchosen vertices are omitted for the display. (d,e) The GBS 
enhancement on finding the 80-vertex dense subgraph out of a 144-vertex graph for SA (d) and RS (e) algorithm. The experiment obtains a 
mean photon-click number of 61. A set of 80-photon-click samples are post-selected for this study. The X-axis corresponds to the number of 
steps used in the optimization process, while the Y-axis indicates the corresponding largest density found. The curve is the maximum (mean) 
largest density out of 120 (20) trials for SA (RS), and the shaded area indicates standard deviation (error). In (d) The horizontal dashed line 
shows the density found by the deterministic greedy algorithm, and it is surpassed by the GBS-enhanced SA algorithm. 


The principle of the GBS enhancement on solving the two experiments we need to consider three realistic derivations. 
problems by stochastic algorithms can be understood from the (i) The sampling matrix retrieved from the experiment is not 
concept of proportional sampling. Since the GBS samples are always ideally non-negative as in the original proposal. Imper- 
more likely to have a larger Hafnian in modulus (hereinafter fection in circuits can introduce negative or imaginary terms 
we use Hafnian to refer to Hafnian in modulus), it also holds into the sampling matrix. (ii) Experimental noise like photon 
that subgraphs corresponding to the GBS samples are more loss causes mixed state sampling in GBS, which can result 
likely to have larger Hafnian. Therefore, one can use the GBS in biased diagonal block matrix A, and nonzero off-diagonal 
samples to boost the effectiveness of stochastic algorithms in block matrix L Æ 0 [46]. (iii) Threshold detectors are usually 
solving the Max-Haf problem by augmenting its success prob- used instead of photon-number-resolving detectors [47]. 


ability. To check whether the proportional sampling mechanism 


Furthermore, it is proved that for a graph of Os and Is, holds for the GBS with threshold detectors on complex-valued 


its density is positively correlated to Hafnian, and dense k- sampling matrices, we perform Monte Carlo simulation to re- 
subgraph problem solving can also be expected to gain en- veal the numerical correlation between Torontonian and Haf- 
hancement from GBS [10]. From another point of view, nian or density for randomly generated complex sampling ma- 
by working in a quantum-classical hybrid scheme, the GBS trices (see the Supplemental Material [26]). Shown in Fig. 
serves as an oracle to significantly narrow down the combina- [i{b) and (c), the positive correlation between Torontonian and 
torial search space of the stochastic algorithm since the sub- Hafnian or density validates the underlying principle of pro- 
graphs of small Hafnian or density are unlikely to be sampled. portional sampling, and portends the occurrence of GBS en- 
These two graph problems differentiate from each other hancement. 
by means of their target function’s computational complex- We proceed to test the GBS enhancement on solving Max- 
ity. Hafnian is hard-to-compute while density can be evalu- Haf problem and dense k-subgraph problem. Two stochas- 


ated efficiently. Investigation on the two graph problems of ġe algorithms—namely, the random search (RS) and sim- 
distinct properties provides us with insights into the depen- ulated annealing (SA) [9] [10J—are studied. RS represents 
dence of GBS enhancement on the computational complexity the naive way of solving the combinatorial problem by uni- 
of the graph feature itself. formly sampling from the whole solution space, which is free 

While the above discussion holds for the ideal GBS, in from being trapped by local optimum but is costly and inef- 
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FIG. 3. The scaling of the GBS enhancement benchmarked by the score advantage and speed advantage as a function of photon-click number. 
The problem is defined on a randomly generated complex-valued full graph of 144 vertices. (a),(b) The GBS enhancement on the Max-Haf 
problem with various photon-click numbers. The mean photon-click number of the experiment is 21. The score advantage as a function of 
photon-click number is shown in (a), which is defined as the ratio of maximum Hafnian in square of modulus searched at 1000 steps by the 
GBS-enhanced RS algorithm to that searched by the RS algorithm. The speed advantage, which is defined as the ratio of the number of steps 
reaching the target value by the RS algorithm to that by the GBS-enhanced RS algorithm, is shown in (b) as a function of photon-click number. 
The target value of each trial is set as that reached by the RS algorithm at 1000 steps. A clear rising trend with increased photon-click number 
can be observed for both the score advantage and the speed advantage. (c),(d) The GBS enhancement on the dense k-subgraph problem for 
various photon-click number. The mean photon-click number of the experiment is 61. The score advantage is displayed in (c), which is 
defined as the ratio of the density optimized at 10000 steps by the GBS-enhanced RS algorithm to that by the RS algorithm, as a function of 
photon-click number. The speed advantage which is defined as the ratio of the number of steps reaching the same density by the RS algorithm 
to that by the GBS-enhanced RS algorithm, is showns in (d) for various photon-click numbers. For each trial, the target value is set as that 
reached by the the RS algorithms at 10000 steps. No significant increasing trending with photon-click number is observed for this problem. 
Error bars indicate standard error. 


ficient. SA combines mechanisms from both random explo- 
ration that prevents it from being stuck in local minima and 
hill climbing that enables it to approach good solutions fast, 
but proper choice of parameters is crucial for guaranteeing the 
algorithm’s performance. Together, the two algorithms of dis- 
tinct working subroutines help benchmark the enhancement 
of GBS on graph applications more comprehensively. 


The experiment is performed on a randomly generated and 
fully connected 144-mode optical interferometer, and a subset 
of samples with coincident photon-click number up to 80 are 
used for the study. Figure [2{a) and (b) show the maximum 
Hafnian of a 12-vertex subgraph on a 144-vertex full graph 
found for the two algorithms and their GBS-enhanced variants 
as a function of searching steps. For both the RS and SA algo- 
rithm, it is evident the GBS-enhanced variants improve the ef- 
fectiveness of the algorithms by finding larger Hafnian within 
the same steps. An illustration of the full graph correspond- 
ing to the experiment, together with the subgraph searched by 
the GBS-enhanced SA algorithm in highlight, is shown in Fig. 


Bto). 
Similarly, Figure BG) and (e) plot the largest density found 


at various steps for the four algorithms. The GBS samples are 
80 photon-click events from the 144-mode quantum device, 
which are in the quantum advantage regime. On average each 
sample would take Frontier, the current fastest supercomputer 
in the world [48], more than 700 seconds to generate using 
exact methods, as estimated with the state-of-the-art classical 
sampling algorithm [49], and we used 221891 samples in total 
for the study which amount to ~ 5 years on Frontier. It is 
observed that both RS and SA algorithms gain enhancement 
from the GBS samples in searching for subgraphs of higher 
density at the given step. Specially, it is noted that the density 
found by the deterministic Greedy algorithm [50], which is 
marked as the horizontal dashed line, can be outperformed by 
the GBS-enhanced SA algorithm, confirming the advantage of 
stochastic algorithms. 


Having established the GBS enhancement, we continue to 
investigate how this enhancement scales on our device. We 
benchmark the GBS enhancement by defining the score ad- 
vantage and speed advantage. The former is that, for a given 
step, the maximal score (in terms of the Hafnian or the den- 
sity) obtained by the GBS-enhanced algorithms divided by 
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FIG. 4. Histogram of steps for GBS-enhanced RS algorithms to achieve the target value as a function of various noise levels. For each subplot, 
a number of trials are repeated, and the y axis corresponds to the counts of the occurrence of the number of steps to surpass the target value. 
The target is set as the mean value reached by the classical RS algorithm at 10000 steps. (a),(b) Theoretical simulation of the effect of photon 
loss noise on the GBS enhancement on the Max-Haf problem (a) and the dense k-subgraph problem (b). The plots show the GBS enhancement 
with samples from a theoretically simulated sampler with no photon loss 7 = 1, and theoretically simulated sampler with system efficiency 
n = 0.75 and 7 = 0.5. It is observed that the ideal sampler without photon loss needs fewer steps for the task than the lossy sampler. (c),(d) 
Simulating thermal noise’s influence on GBS-enhancement for (c) the Max-Haf problem and (d) the dense k-subgraph problem. The plots 
exhibit the GBS enhancement with samples from a theoretically simulated sampler of zero thermal noise e = 0, and theoretically simulated 
samplers of thermal noise level e = 0.25 and e = 0.5. The results indicate that the simulated sampler with no thermal noise needs fewer steps 
to achieve the target value. (e),(f) The GBS-enhancement with samples from a low noise level experiment (red) and a controlled high noise 
level experiment (brown) on (e) the Max-Haf problem and (f) the dense k-subgraph problem. Fewer steps are needed by the low noise level 
experiment than the high noise level experiment. The low noise level experiment corresponds to an averaged system efficiency 7 = 0.472 and 
averaged thermal noise level € = 0.02, whereas the high noise level corresponds to an averaged system efficiency 7 = 0.333 and an averaged 
thermal noise level € = 0.192. For all plots the problem is defined on a 144-vertex full graph with a photon click of 10. The low- (high-) 
noise level experiment has a mean photon number of 7 (6), and samples of photon click 10 are used for the investigation. Error bars indicate 
statistical fluctuation. 


that by the classical algorithms only. The latter is that, to reach no obvious increasing trend at larger size. A similar be- 


a target score, the ratio of the needed searching steps between 
the GBS-enhanced and classical algorithms [26]. We use the 
parameter-free RS algorithm to probe the scaling properties. 

Figure Bla shows the scaling of the score advantage of the 
GBS enhancement for a fixed 10° steps. Remarkably, the 
score advantage rises up steadily, from ~ 24 at a photon click 
of 12 to ~ 92 at a photon click of 28. The speed advantage is 
plotted in Fig. [3{b) also as a function of photon-click number. 
The speed advantage starts at ~ 89 at a photon click of 12 
and becomes increasingly larger as the size increases, reach- 
ing ~ 212 at 28 photons. Here, due to the computational over- 
head, we use the Sunway TaihuLight supercomputer to eval- 
uate the Hafnian. Overall, the results of Fig. [3{a,b) provide 
strong evidence that the GBS enhancement as benchmarked 
by the score advantage and speed advantage increases with 
the photon-click number in solving the Max-Haf problem by 
RS algorithm on Jitizhdang. 

Figure Bto) plots the score advantage for the dense k- 
subgraph problem at increasing photon-click number. While 
all the data points show a positive advantage (>1), there is 


havior is observed in the speed advantage, as shown in Fig. 
[3{d). In the Supplemental Material we show numerically 
simulated results for an ideal GBS sampler on the dense k- 
subgraph problem of both randomly generated non-negative- 
valued graph and complex-valued graph, which exhibit trend- 
ing of the score advantage and speed advantage similar to that 
reported in our experiment. 


Noise is a major problem for the NISQ device. In the 
GBS, photon loss (which can be caused by limited ef- 
ficiency of the optical elements and detection) and thermal 
noise (which can be caused by spatial mode mismatch 
of the sources) can turn the pure-state GBS into mixed-state 
GBS. For the graph problem solving, these noises can make 
the sampling matrix deviate from the ideal [29], and could 
decrease the positive correlation between Hafnian or density 
of the encoded matrix and that of the sampling matrix. To 
characterize the influence of these noises on GBS enhance- 
ment on the graph problem solving, we benchmark them with 
the RS algorithm that is free from parameter choosing. We 
compare the steps needed for achieving a target value of the 


problem between samplers of various noise levels. The prob- 
ability distribution of the steps follows the Geometric Distri- 
bution, which gives the probability that the first occurrence of 
success requires k independent trials. 


P(X =k) =(1—p)*""p (5) 


where k = 1,2... is the number of steps, and p is the proba- 
bility that a GBS sample can produce a better result than the 
target. The noise’s influence on the GBS enhancement can 
be simply benchmarked by the parameter p, since a larger p 
would indicate fewer steps are needed which correspond to a 
stronger GBS enhancement, and vice versa. 

To investigate the effect of photon loss, we theoretically 
simulate the performance with an ideal sampler and that with 
an overall photon loss of 25% and 50%, for the same optimiza- 
tion task. Figure [4] (a) and (b) show histograms of the num- 
ber of steps for the GBS-enhanced RS algorithms to achieve 
the target value. There is a significant reduction of the steps 
at increasing system efficiency 7. The p value of the sam- 
pler with a unit efficiency is 0.0196 (0.0024) for the Max-Haf 
(dense k-subgraph) problem, whereas the lossy sampler with 
efficiency 7 = 0.5, 0.75 correspond to a p = 0.0039 (0.0012), 
p = 0.0071 (0.0018), respectively. The results indicate that 
lower photon loss will lead to a stronger GBS enhancement. 
A recent theoretical study reported similar findings, including 
the dependence of the GBS enhancement on the partial photon 
distinguishability [51]. 

Figure[4](c) and (d) show the theoretical simulation results 
for the thermal noise [52]. Three examples are studied, where 
the thermal noise is chosen for €e = 0,0.25,0.5. Again, a 
strong decrease of the required steps is observed for lower 
thermal noise. The p value of the ideal sampler for the Max- 
Haf (dense k-subgraph) problem is 0.0194 (0.0024), whereas 
the p value of the sampler with 0.25,0.5 thermal noise is 
0.0023 (0.0011), and 0.0008 (0.0006), respectively. The re- 
sults show the importance of eliminating the thermal noise to 
achieve a higher GBS enhancement. 

Having studied these effects theoretically, we now bench- 
mark the noise influence experimentally. The experiment re- 
sults at a typical noise level, where 7 = 0.472 and € = 0.02, 
is compared with a controlled higher noise level, 7 = 0.333 
and € = 0.192. As shown in Fig. [4fe) and (f), the exper- 
imental samples with low noise level demonstrate stronger 
GBS enhancement for both graph problems, which is in good 
agreement with the theoretical simulation. The p value for 
the low noise experimental sampler on the Max-Haf (dense 
k-subgraph) problem is 0.0031 (0.0012), whereas it is 0.0014 
(0.0008) for the controlled high noise experimental sampler. 
Interestingly, samples from the modest-noise-level experi- 
ments, though with less enhancement, can still improve RS 
algorithm. 

In this Letter, we have demonstrated the GBS enhancement 
on stochastic algorithms in solving two graph problems of 
distinct properties with the 144-mode NISQ device Jilizhdng 
in the quantum computational advantage regime. It is an 


open question, however, whether the GBS can yield advan- 
tage compared to improved classical algorithms and quantum- 
inspired algorithm. Also, the GBS enhancement can depend 
on the properties of the input graphs, for which more com- 
prehensive algorithm analysis and discussions for various sit- 
uations are expected. We hope that our work will stimu- 
late experimental efforts on larger-scale, higher-fidelity and 
fully programmable GBS, exploration of real-world applica- 
tions where the computational problems can be mapped onto 
the GBS, and development of more efficient classical and 
quantum-inspired algorithms. 
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